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Abstract 

In this article, we present a family of numerical approaches to 
solve high-dimensional linear non-symmetric problems. The principle 
of these methods is to approximate a function which depends on a large 
number of variates by a sum of tensor product functions, each term of 
which is iteratively computed via a greedy algorithm |20] . There exists 
a good theoretical framework for these methods in the case of (linear 
and nonlinear) symmetric elliptic problems. However, the convergence 
results are not valid any more as soon as the problems considered are 
not symmetric. We present here a review of the main algorithms pro- 
posed in the literature to circumvent this difficulty, together with some 
new approaches. The theoretical convergence results and the practical 
implementation of these algorithms are discussed. Their behaviors are 
illustrated through some numerical examples. 

Introduction 

High-dimensional problems arise in a wide range of fields such as quantum 
chemistry, molecular dynamics, uncertainty quantification, polymeric fluids, 
finance... In all these contexts, one wishes to approximate a function u 
depending on d variates x\, where d G N* is typically very large. 

Classically, the function u is denned as the solution of a Partial Differen- 
tial Equation (PDE) and cannots be obtained by standard approximation 
techniques such as Galerkin methods for instance. Indeed, let us consider a 
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discretization basis with N degrees of freedom for each variate (N G N*), so 
that the discretization space is given by 

V N := Span {^Oi) ■ • • ^f(x d ), 1 < i 1} ■ ■ ■ , i d < iv} , 
where for all 1 < j < d, (ip^ M is a family of A" functions which only 

\ / l<i<N 

depend on the variate Xj. A Galerkin method consists in representing the 
solution u of the initial PDE as 

u(x u ■■ - ,x d ) » X h,- M^h ( x i) ■ ■■4f( x d), 

l<ii,- ,id<N 

and computing the set of N d real numbers (Aj lr .. ,i d ) 1<il ... i d <Ar- Thus, the 
size of the finite-dimensional problem to solve grows exponentially with the 
number of variates involved in the problem. Such methods cannot be imple- 
mented when d is too large: this is the so-called curse of dimensionality [2]. 

Several approaches have recently been proposed in order to circumvent 
this significant difficulty. Let us mention among others sparse grids |21) . ten- 
sor formats jll] , reduced bases j3] and adaptive polynomial approximations 

i- 

In this paper, we will focus on a particular kind of methods, originally 
introduced by Ladeveze et al. to do time-space variable separation [12], 
Chinesta et al. to solve high-dimensional Fokker-Planck equations in the con- 
text of kinetic models for polymers [TJ and Nouy in the context of uncertainty 
quantification |15j . under the name of Progressive Generalized Decomposition 
(PGD) methods. 

Let us assume that each variate Xj belongs to a subset Xj of M. mj , where 
rrij E N* for all 1 < j < d. For each d-uplet (r^,-- - , r^) of functions 
such that r^ only depends on Xj for all 1 < j < d, we call a tensor product 
function and denote by r^ 1 ' <g> • • • r^ d ' the function which depends on all the 
variates x±, ■ • ■ ,x d and is defined by 

m (A\ X-\ x ■ ■ • x X d — y IR 

{ (x 1 ,---,x d ) i y r w {xx) ■ ■ -r w (x d ). 

The approach of Ladeveze, Chinesta, Nouy and coauthors consists in 
approximating the function u by a separate variable decomposition, i.e. 

n n 

u(x ir - • ,x d ) « ^Vj^Oi) • • -r^Zd) = • -(g)rj, d) (a;i, • • • ,x d ), (1) 

fc=i fe=i 
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for some n G N*. In the above sum, each term is a tensor product function. 
Each d-uplet of functions \ r^\ ■ ■ ■ ,rj®J is iteratively computed in a greedy 

|2U] way: once the first k terms in the sum (jTJ have been computed, they 
are fixed, and the (k + l) th term is obtained as the next best tensor product 
function to approximate the solution. This will be made precise below. 

Thus, the algorithm consists in solving several low-dimensional problems 
whose dimensions scale linearly with the number of variates and may be 
implementable when classical methods are not. In this case, if we use a 
discretization basis with TV degrees of freedom per variate as above, the 
size of the discretized problems involved in the computation of a d-uplet 
(r^\ • • • , rjf^j scales like Nd and the total size of the discretization problems 
is nNd. 

This numerical strategy has been extensively studied for the resolution of 
(linear or nonlinear) elliptic problems (2"U| [T31 EH El HE]- More precisely, let 
u be defined as the unique solution of a minimization problem of the form 

u — argmin£(/z;), (2) 
vev 

where V is a reflexive Banach space of functions depending on the d variates 
Xi, Xdi and £ : V — > R is a coercive real- valued energy functional. Besides, 
for all 1 < j < d, let V Xj be a reflexive Banach space of functions which only 
depend on the variate Xj. The standard greedy algorithm reads: 

1. set uq = and n — 1; 

2. find (rn \ ■ ■ ■ , r$\ G V X1 x • ■ ■ x V Xd such that 

{ r n\ ■■■ i r n ] ) e argmin £ (u n _i + r (1) ® • • • ® r (d) ) , 

(rW.-.rWjeF^x-xK,, 

3. set n n = n n _i + r^ 1 ^ (gi • ■ ■ (g) and n = n + 1. 

Under some natural assumptions on the spaces V, V^ i; T4 d and the energy 
functional £, all the iterations of the greedy algorithm are well-defined and 
the sequence (u n ) n6 N* strongly converges in V towards the solution u of the 
original minimization problem (j2J). 

This result holds in particular when u is defined as the unique solution of 

{find u G Usuch that 
Vv G V, a(u,v) = l(v), 
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where V is a Hilbert space, a a symmetric continuous coercive bilinear form 
on V x V and / a continuous linear form on V. In this case, u is equivalently 
solution of a minimization problem of the form ([2]) with £(v) = \a{v, v)—l{v) 
for all v G V. 

However, when the function u cannot be defined as the solution of a 
minimization problem of the form designing efficient iterative algorithms 
is not an obvious task. This situation occurs typically when u is defined as 
the solution of a non-symmetric linear problem 

{find u G Vsuch that 
Vi> G V, a(u,v) = l(v), 

where a is a non-symmetric continuous bilinear form oiil^x^ and I is a 
continuous linear form on V. 

The aim of this article is to give an overview of the state of the art of 
the numerical methods based on the greedy iterative approach used in this 
non-symmetric linear context and of the remaining open questions concern- 
ing this issue. In Section [TJ we present the standard greedy algorithm for 
the resolution of symmetric coercive high-dimensional problems and the the- 
oretical convergence results proved in this setting. Section [2] explains why a 
naive transposition of this algorithm for non-symmetric problems is doomed 
to failure and motivates the need for more subtle approaches. Section [3] de- 
scribes the certified algorithms existing in the literature for non-symmetric 
problems. All of them consist in symmetrizing the original non-symmetric 
problem by minimizing the residual of the equation in a well-chosen norm. 
However, depending on the choice of the norm, either the conditioning of the 
discretized problems may behave badly or several intermediate problems may 
have to be solved online, which leads to a significant increase of simulation 
times and memory needs compared to the original algorithm in a symmetric 
linear coercive case. So far, there are no methods avoiding these two prob- 
lems and for which there are theoretical convergence results in the general 
case. In Section HI we present some existing algorithms designed by Nouy [16] 
and Lozinski [T4] to circumvent these difficulties and the partial theoretical 
results which are known for these algorithms. Section is concerned with 
another algorithm we propose, for which some partial convergence results 
are proved. In Section [61 the behaviors of the different algorithms presented 
here are illustrated on simple toy numerical examples. Lastly, we present in 
the Appendix some possible tracks to design other methods, but for which 
further work is needed. 
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1 The symmetric coercive case 



1.1 Notation 

Let us first introduce some notation. Let d be a positive integer, mi, md 
positive integers and X\, ...,X d open subsets of IR™ 1 , IR md respectively. 
Let fi Xl , ■■■■> l^x d denote measures on X\, Xd respectively. Let 1?{X\, fi Xl ), 
L 2 (Xd] /J> Xd ) be associated L 2 spaces, i.e. vectorial spaces which are com- 
plete when endowed with the scalar products 

V/,0 G L 2 {Xj] fj, Xj ), (f,g) Xj := / f{x j )g{x j ) /jL x .(dxj), VI < j < d, 

Jx, 



and their associated norms || ■ || • \\x d - For instance, in the case when 

X\ = (0, 1) and /j, Xl is the standard Lebesgue measure on Xi, the spaces 

L 2 (0, 1), L 2 pcr (0, 1) and L 2 (0, 1) := {/ e L 2 (0, 1), f = o} are examples of 

such L 2 spaces. 

In the rest of this article, for the sake of simplicity, we will omit the 
reference to the measures /i xi , fj, Xd and denote by L 2 (Xi) = L 2 (Xi] fJ, xl ), 
L 2 (X d ) = L 2 (X d ; fi Xd ). 

We introduce the space L 2 {X X x • • • x X d ) := L 2 {X X ) ® ■ ■ ■ ® L 2 (X d ). This 
space is a Hilbert space when endowed with the natural scalar product 

Vf,g G L 2 {X x x- ■ -xX d ), (f,g}:= / f(x x , ■ ■ ■ , x d )g{x 1} ■ ■ ■ , x d ) n xi {dxi) ■ ■ ■ n Xd {dx d ), 

JXix-xX d 

and the associated norm is denoted by || • H^x-x^- 

Let V C L 2 (X 1 x • • ■ x Xd), V X1 C L 2 {X X ), V Xd C L 2 (X d ) be Hilbert 
spaces endowed respectively with scalar products denoted by (•, -)y, (-, -)v xi i 

(-,-}v Xd and associated norms || • || y , || ■ \\ Vxi , || ■ \\v Xd - 

We define V, V' xv V' Xd as the dual spaces of V, V X1 , V Xd with 
respect to the L 2 scalar products (-, ■), (•, ■)x 1 , (■, -)x d - These dual spaces 
are endowed with their natural norms || • ||y/ etc. 

Lastly, the Riesz operator Ry : V — > V is defined by 

\/v,w G V, (v,w) v = (R v v,w)v,v- 
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It holds in particular that ||i>||v = Similar operators Ry xi , •••> Rv x 

are introduced for the spaces V X1 , V Xd . 



For any <i-uplet (r^ 1 ', ■ • • , r^) G V Xl x • • • x V Xd , we define the tensor 
product function r^ 1 ' (g) • • • (g) as follows 



ni mi <*vi x • • • x Afw — ^ ffi 

r (1) ® ■ • • ® r (d) : <{ / A d au . fdl , . 



In the particular case when d = 2, we shall denote respectively Xi, A4, 
mi, by x, X, m x , V x and x 2 , X 2 , m 2 , 14 2 by t, T, m t , V t . 

Besides, for any Banach spaces Hi, H 2 , the space of bounded linear op- 
erators from Hi to H 2 will be denoted by £(Hi, H 2 ). 

1.2 Theoretical results 

We recall here the theoretical framework of the standard greedy algorithm 
in the coercive symmetric case. 

Let us consider the problem 



where 

• a(-, •) is a symmetric, coercive continuous bilinear form on V x V; 

• I is a continuous linear form on V. 

Then, problem ([3]) is equivalent to the minimization problem 




(3) 




(4) 



where 





The greedy algorithm reads: 



1. let Uo = and n = 1; 
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2. define \Tn\ • • • , E V X1 X • • • x V Xd such that 
( r n 5 " ' j r i d) ) G argmin £ + 

(r( 1 ),-,rW)eV a!1 x-xV a!d 



(6) 



3. define u n = ii n _i + <E> • • • ® and set n = n + 1. 
Let us denote by 

£ := { r W®--.®rW, r« G , • • ■ G V Xd } (7) 
and make the following assumptions: 



(Al) Span(S) = V; 

(A2) S is weakly closed in V. 

These assumptions are usually satisfied in the case of classical Sobolev spaces [5]. 

Theorem 1.1. Assume that (Al) and (A2) are satisfied. Then, for all n E 
N*, there exists at least one solution \ Tn \ • • • , r^^\ E V Xl x • • • x V Xd (not 

necessarily unique) to fifty and any solution satisfies Tn ®- • -®Tn ^ if and 
only i/ «„_i 7^ u. Besides, the sequence (w n )neN* strongly converges towards 
u in V . 

The following Lemma will be used later. Although the proof is given 
in 1101. we recall it here for the sake of self-containedness. 



Lemma 1.1. For all v E V, let us denote by \\v\\ a := ^Ja(v,v). Then, for 
all n EW, 



® ■ ■ ■ ® r^W = 



a(u — m„_i, (g) • • • <g) r^] 



(rCi),...,rW)eT4 1 x-x\4 d ,rCi)®-®rW^0 ll r(1) ® " " " ® r(d) l 

(8) 

Proof. Let us prove ([H]) for n = 1 . The proof is similar for larger n£N*. The 
(i-tuple ^r[ X \ • • • , ~\ E V X1 x • • • x V^. d solution of (|6]) for n = 1 equivalently 
satisfies: 

r^\ ■ ■ ■ , G argmin - llu — (g) ■ ■ • ® || 2 . (9) 

y (r(i) r , r M) 6 v lI x...x^ 2 



The Euler equations associated to this minimization problem read: for all 
(5 r M ... >M) eVxi x •••xK, 



a frf 5 (8) • • • <g } , rf } <g • ■ ■ <g rf* X) <g 5r {d) + rf } <g • • ■ <g rf* 2) <g 5r (d_1) <g rf H h 5r (1) (g 

r{ 1} <g • • • g) rj d_1) ® <5r w + r[ 1] <g • • • ® rf~ 2) <g tfr^" 1 ) <g + • • • + 5r (1) g> rf } <g • • • <g 



a ii 



which implies that 



„tf) 



a l w, 



Let now (r^\ • • ■ , r^) G T4 1 x • • ■ x V Xd be such that g) 
Using (ED and (HO]), it holds that 



(10) 
r (d) ^ o. 



a I u, (g ■ ■ • <g ?i 



u 



.0) 



<g • • • <g ?i 



2 ' i 



u — r 



(i) 



„(<*) 



< 



(u, g) • • • (g r ( 
||r( x ) (g) • • • <gi 



Therefore, 



\ 2 

J a (u, 



(g • ■ ■ (g 



|rW (g ■ • • <g r^\\ a 



Taking the supremum over all [r^\--- , r^) G x ••• x V Xd such that 
(g • • ■ (g ^ yields the result. □ 



Equation (jSJ) implies in particular that for all n G N* 



(g • • ■ r, 



(d)\ 



sup 

r m,-,r( d ))eV Xl x-xV x 



r (d) ) - a (u n -i, r (1) g) • • ■ <g r (d) ) 



|rW <g • ■ • (g r( d ) | 
(11) 



Let us rewrite the greedy algorithm in the particular case when d = 2. 
1. Let Uq = and n — 1; 
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2. define (r n , s n ) G V x x V t such that 

(r n ,s n )e argmin £ (w n _i + r <g) s) ; (12) 

3. define ?/ n = u n -\ + r n ® s n and set n = n + 1 . 

For the sake of simplicity, in the rest of the article, all the algorithms will 
be presented in the case when d = 2. The generalization of the approaches 
to a larger number of variates d is straightforward unless mentioned. 

The Euler equations associated to the minimization problem ( JT2|) read 

a(u n -x + r n ®s n ,5r<g!S n + r n <$8s) = l(5r®s n + r n ®5s), V(5r,5s) G 14 x V t . 

(13) 

As a consequence of Theorem II. 1[ provided that the set 

S = {r ® s, r G 14, s G V*} (14) 

satisfies assumptions (Al) and (A2), at the first iteration of the algorithm 
(n — 1), as soon as the form I is nonzero, there exists at least one solution 
(ri.si) G 14 X V t of 

a(r\ (g> si, 5r <g> si + T\ <E> 5s) = /(5r (g) si + n <8> 5s), V(5r, 5s) G 14 x Vt, 

such that ri ® Si 7^ 0. 

In practice, at each iteration n G N*, a pair (r n , s n ) G 14 X Vt is computed 
via the resolution of the Euler equations (TLB"]) using a fixed-point procedure 
which reads as follows: 

• choose (^Tn^ , Sn^ ) G 14 x V t and set m = 1 ; 

• find (Vi m) , s^ m) ) G 14 x 14 such that 

a (u n ^ + ri m) ® s<T 1} , 5r ® si" 1 " 1 )) = Z (<5r ® s^T 1} ) , V5r G 14, 
a + ri m) ® si m) , ri m) ® 5s) = / (ri m) ® 5s) , V5s G V t ; 

(15) 

• set m = m + 1. 

This fixed-point algorithm is numerically observed to converge exponen- 
tially fast in most situations, although, at least to our knowledge, there is no 
rigorous proof in the general case. 
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2 The non- symmetric case 
2.1 General framework 

Let us now consider the case of a non-symmetric linear problem of the form 

{find ueV such that , . 

Vv G V, a(u,v) = l(y), ^ ' 

where 

• a(-, ■) is a nonsymmmetric continuous bilinear form on V x V; 

• I is a continuous linear form on V. 

In the rest of the article, we will assume that 

(A3) problem ([IE]) has a unique solution u G V for any continuous linear 
form / G £(V,R). 

We denote by .A G £(V, V) the operator defined by 

Vi>, w E V, (Av, w)y = a(v, w), 

and by L the element of V such that 

Vi;GF, (£,v)v = Z(u). 

We also introduce the operator A : V — > V and the linear form L EV 1 
defined by A = RyA and L = RyL so that the unique solution u to (fT6|) is 
also the unique solution to the problem 

{find u G V such that 
Au = L in V. 

It follows from assumption (A3) that A and A are invertible operators. 
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2.2 Prototypical examples 

Let us present two prototypical examples we will refer to throughout the rest 
of the paper. 

• The first one is 

f find u G Hl(X) ® L 2 (T) such that , , 

\ -A x u + b x ■ W x u + u = f in V'{X x T), ^ ' 

with / G if -1 ^) <g> L 2 (T) and 6 X G M m *. For this problem, V = 
Hq(X) ® L 2 (T), V" = tf-V) ® L2 C7") and 

Vw, w G V, a(w, = f XxT (y x u ■ V x v + v(b x ■ V x u) + uv) , 
Wv G V, l(v) = J T (f,v) H -i (x)tH i {x) . 

In this case, A = —A x + b x ■ V x + 1. 

• The second example is 

find u e Hq(X x T) such that , v 

-A Xit u + 6 • V^u + u = f in x 71, ( j 

with / G H~ 1 (X x T) and 6 = & t ) G M. mx x M mt . For this problem, 
V = Hl{X xT),V' = H"\X, T) and 

Wu, v G V, a(u, v) = f XxT (V x ju ■ V Xjt v + v(b ■ V^it) + uv) , 
Vt; G V, l(v) = (f,v) H -i (XxT)>H i {XxT) . 

In this case, A = —A xt + b ■ V Xj t + 1. 



2.3 Failure of the standard greedy algorithm 

Problem (ITtjj) cannot be written as a minimization problem of the form (j3J) 
with an energy functional given by ([5]) . The definition of the greedy algorithm 
via the minimization problems Q) or f|T2|) cannot therefore be transposed 
to this case. However, a natural way to define the iterations of a greedy 
algorithm for the non-symmetric problem 0161) is to define iteratively for 
n G N* the pair (r n , s n ) G V x x V t as a solution of the following equation 

a(u n _i +r„ eg) s n , <5r <g> s n + r n <g> 5s) = /(<5r <g> s n + r n CgxJs), V(5r, 5s) G V x xV t , 

(19) 
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by analogy with the Euler equations (II 3 p . This is the so-called PGD-Galerkin 
algorithm [3J. 

Actually, there are cases when I ^ and any solution (r±, si) G V x x V t 
of the first iteration of the algorithm 

a(n <g> si, 5r <g> si + n <g> 5s) = l(8r <g> si + r x <g> 5s), V(£r, 5s) G K x V t , (20) 

necessarily satisfies n <g) Si = 0. Such an algorithm cannot converge since 
the approximation u n = (g) s^ given by the algorithm is equal to for 

k=l 

any n G N*. Besides, this situation may occur even when the norm of the 
antisymmetric part of the bilinear form a(-, ■) is arbitrarily small. 

Let us give an explicit example. 

Example 2.1. Let X — J' — (—1,1) and /j, x (respectively /j, t ) be the Lebesgue 
measure on X (respectively on T). Let b G R, V x = iZj^-l.l), V t = 
L 2 (— 1, 1) and V = V x ®V t . Consider the non-symmetric problem l[Tb)) with 



Wv, w G V, a(v, w) 



(V x v ■ V x w + (6 ■ V x v)w + vw) 



XxT 



and 



Wv G V, l(v) = [ fv, 

J XxT 



with f G L 2 peT {-l,l) ® L 2 (-l,l). 
Problem fT(j\) is equivalent to 



find u G H^ er (-1, 1) <g> L 2 (-l, 1) such that 
-A x u + bV x u + u = f mD'(lxT). 



(21) 



In this context, equations §Wy read 

find (ri, si) G ifp er (— 1, 1) x L 2 (— 1, 1) such that 
/i ki(t)| 2 rft] (-r?(x) +K(x) + n(x)) = J' 1 f(x,t)s 1 (t)dt, 
f^iKWl 2 + \ r i(x)\ 2 ) dx Si(t) = f^ 1 f(x,t)r 1 (x)dx, 



(22) 



smce £/ie periodic boundary conditions on r\ imply that J_ x r 1 (x)r' 1 (x) = 0. 
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Unlike the symmetric case, there exists an infinite set of functions f G 
Lp er (— 1, 1) Cg> L 2 {— 1, 1) such that f ^ and any solution (r l5 s\) £ V x x V t 
of equations |HP necessarily satisfies r\ (8> S\ — /or any arbitrarily small 
value of \b\. This is the case for example when f(x,t) = <p(x — t) for all 
(x,t) elx (—1, 1) with (p G Lp Cr (— 1, 1) an odd real-valued function. 

Let us argue by contradiction. If (ri,si) G V x x V t is a solution to (2B) 
such that n <E> Si 7^ 0, up to some rescaling, we can assume that 

J \ Sl (t)\ 2 dt = j 1 (\r[(x)\ 2 + \n(x)\ 2 ) dx = A > 0. 

Thus, we can rewrite |Hj) as 

—r"(x) + br[(x) + r±(x) 

Sl (t) 



f(x,t) Sl (t)dt, 



f(x, t)r 1 (x) dx. 



Plugging the second equation into the first one, we obtain 

•1 / pi 



- r'l(x) + br[(x) + n(x) = J (J 



f(x,t)f(y,t)dt)r 1 (y)dy. (23) 



Let us denote by g(x, y) = f(x, t)f(y, t) dt for all (x, y) G IR 2 . As (j) is an 
odd, 2-periodic function, it holds that 

g{x,y) = J f(x,t)f(y,t)dt 

(j)(x - t)<j)(y - t) dt 
-1 

<p(x - t)(f>(t - y) dt 

-1 

<p(x — y — u)<p(u) du 

-l+y 

•1 

4>(x — y — u)<p(u) du. 

-1 
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Taking the Fourier transform of equation < fl?5]) yields that for all k G 7rZ 7 
(|fc|2 + lhk + 1)fl(A;) = _i_ Vi(A;), 

u>aere 

1 Z" 1 

— 2 y ?"i(a;)e lk ' x dx. 
Futhermore, A G and 0(0) = an odd function). Thus, since 



(f>(k) is a purely imaginary number, — yp(k)j = 4>{k) and a solution n 
necessarily satisfies r\{k) = for all k G 7rZ, which yields a contradiction. 

This example clearly shows that a naive transposition of the greedy algo- 
rithm to the non-symmetric case by analogy with the Euler equations ( Jl3|) 
obtained in the symmetric case may be doomed to failure. 

This article presents a review of some methods which aim at circumvent- 
ing this difficulty A particular highlight is set on the practical implementa- 
tion of these methods and on the existence of theoretical rigorous convergence 
results. The properties of the different algorithms which are dealt with in 
this article are summarized in Figure 1. 



3 Residual minimization algorithms 

In this section, we present some numerical methods used for the computation 
of separate variable representations of the solution of non-symmetric prob- 
lems, for which there are rigorous convergence proofs. A natural idea is to 
symmetrize (fT6|) using a reformulation as a residual minimization problem 
in a well-chosen norm. These algorithms are also called Minimum Residual 
PGDin the literature [3]. 



3.1 Minimization of the residual in the L?(X x T) norm 

Let us assume that L G L 2 (X x T) and that there exists D(A) C V a dense 
subdomain of L?{X x T) such that A(D(A)) C L?{X x T). The mapping 
A : D(A) — y L 2 (X x T) defines a linear operator on L 2 (X x T). Let us 
assume moreover that A is a closed operator. This implies in particular that 
D(A), endowed with the scalar product 

Wv,w G D(A), (v,w) D{A ) = (v,w) + (Av,Aw), 
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THEORY 


PRACTICE 


L 2 residual minimization 


Additional regularity on the right-hand side 
is needed though. 


The conditioning of the resulting problems 
scale quadratically with the conditioning 
of the original problem. 


Dual norm residual 
minimization 


OK 


Need to solve several small— or high— dimensional 
symmetric coercive problems 


PGD-Galerkin 


There are cases where the algorithm does not 
converge towards the true solution. 


OK 


Minimax 


OK for separated operators 


OK 


Dual Greedy 


OK in finite dimension when V = V x ® Vt 
Problems in infinite dimension 


OK 


X-Greedy 


Same situation as the Dual-Greedy 


Not clear how to implement the algorithm 
in practice. 


Decomposition 


OK in finite dimension provided that the 
explicit part of the bilinear form is small enough 
compared to its implicit part. 


OK but slow 
Diverges if the explicit part of the bilinear 
form is too large. 



Figure 1: Summary of the different greedy algorithms used for non-symmetric 
high-dimensional linear problems. 
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is a Hilbert space. 

A first approach, inspired by [9], consists in applying a standard greedy 
algorithm on the energy functional 

£(v) = \\Av - L\\ 2 L2{Xxr) , Vv E D(A). 
Let us consider the case when 

A = Y j Af®Af ) 

i=i 

where for all 1 < i < p, Ax and A^ are operators on L 2 (X) and L 2 (T) 
with domains D (^Ax^ and D ^A^^ respectively. We denote by D x = 

f)i =1 D^Ax^ and D t = {~)? =1 D (/^tj > an d assume that D a . and D t are 

dense subspaces of L 2 (X) and L 2 (T) respectively and are Hilbert spaces, 
when endowed with the scalar products 



p 

Vv, «; G A,, (v, w) Da: = (v, w) x + ]T (4V A®w) x 

i=l 

and 

V«, «; G A, («, ^)d 4 = (v, w>r + X] v > A * 

The greedy algorithm reads: 

1. let u = and set n = 1; 

2. define (r n , s n ) G A x A such that 



r 



n ,s n )G argmin ||^4(u n -i + r <g) s) - L\\ 2 L 2 {XxT) ] (24) 

(r,s)eD :c xD t 

3. set u n = w n -i + r n ® s n and n = n + 1. 

Let us denote by S D := {r <g) s, r G A, s G A}- From Theorem ll.l[ 
provided that 

D(A) 



(Bl) SpanS D v ' = D(A); 



1(3 



(B2) S D is weakly closed in D(A); 

the sequence (w n )neN* strongly converges towards u in D(A). 

In the case of problem (fT7|) . A = A x ® A t with A x = —A x + b ■ \7 X + 1 and 
A t = 1, D(A) = (H 2 (X)nHZ(X))®L 2 (T), D x = D(A X ) = H 2 (X) n Hq(X) 
and A = D(A t ) = L\T). 

For problem (US}, A = y^® A^+A^® aJ 2 ) with = -A a .+6 a .-V a .+l, 
A\ l) = 1, Ai 2) = 1 and Af ] = -A t + b t -V t , D(A) = H 2 {X x T)nH^(X x T), 
D x = i/ 2 (A') n tfo 1 ^) and A = H 2 {T) n f^(T). 

In both cases, assumptions (Al) and (A2) are satisfied. 

Actually, when L is regular enough, i.e. if L G ,0(^4*), where A* de- 
notes the adjoint of A and D(A*) its domain, this method is equivalent to 
performing a standard greedy algorithm on the symmetric coercive problem 

A*Au = A*L. 

The Euler equations associated to the minimization problems (|24|) read 
(A(u n ^ + r n ® s n ) - L, A(Sr ®s n + r n (g) 5s)} = 0, V(<Jr, 5s) e D x x D t . 

This method suffers from several drawbacks though. Firstly, the right- 
hand side L needs more regularity than necessary for problem (1161) to be 
well-posed (we need L e L 2 (X x T) instead of L e V'). 

Secondly, and more importantly, the conditioning of the associated dis- 
cretized problems behaves badly since it scales quadratically with the condi- 
tioning of the original problem Au = L. 

3.2 Minimization of the residual in the dual norm 

In order to avoid the conditioning problems encountered when minimizing 
the residual in the L 2 (X x T) norm, another method consists in performing 
a greedy algorithm on the energy functional 

£{v) = \\Av - L\\ 2 V , = H^y 1 ^ - L)\\ 2 Vl \/v G V. 

Here, the residual Av — L is evaluated in the dual norm || • \\y>. In this method, 
the right-hand side L does not need to be more regular than L G V and this 
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approach is equivalent to performing the standard greedy algorithm on the 
symmetric coercive problem 

A*{Ry)~ x Au = A*(R v y 1 L. 

The conditioning of the resulting problem scales linearly with the condition- 
ing of the original Au = L problem. 
The algorithm reads: 

1. let uq = and n — 1; 

2. let (r n , s n ) G V x x V t such that 

(r n , s n ) e argmin \\{R v )~ l [A{u n ^ + r®s)-L] \\ v ; (25) 

(r,s)eV x xV t 

3. set u n = u n ^i + r n ® s n and n = n + 1. 

Provided that E defined by (THj) satisfies assumptions (Al) and (A2), we 
infer from Theorem 1 1 . 1 1 that the sequence (u n ) ne ^ strongly converges to u in 
V. 

The Euler equations associated with the minimization problems fT25|) read: 
for all (8r,8s) G V x x V*, 

(Ry 1 [A(w n _i + r n g) s„) - L] , .Ry 1 [A(Sr ® s n + r„ ® <5s)]) y = 0, 

or equivalently, 

<A(w„_ 1 + r n ® s„) - L, i?^ 1 [A(5r ® s n + r n ® <5s)]) y , y = 0. 

However, even if the conditioning problem of the previous method is 
avoided, this algorithm still requires the inversion of the operator Ry. 

In the case when V = V x ® Vt, the dual space V satisfies V = V' x ® V[, 
so that the operator Ry = Ry x £g> Ry t is a tensorized operator and Ry 1 = 
Ry 1 <S> Ry ■ A prototypical example of this situation is given in problem 
(H2), where we have V' x = H-\X), V( = L 2 (T), R Vx = -A x , i?y t = 1 
and Ry = Ry x ® i?y t . Thus, i?^ 1 = (— A x ) _1 ® 1 and carrying out the 
above greedy algorithm requires the computation of several low- dimensional 
Poisson problems, which remains doable but increases the time and memory 
needs compared to a standard greedy algorithm in the symmetric coervive 
case where b = b x = 0. 
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The situation is even more intricate when V V x ® V t , since the operator 
Ry is not a tensorized operator in general. A prototypical example of this 
situation is problem (|18|) where V = H~ 1 (X x T), Ry = —A x>t and R v l 
cannot be expanded as a finite sum of tensorized operators. These inter- 
mediate symmetric coercive high-dimensional can be solved with a standard 
greedy algorithm presented in Section [T], but this considerably increases the 
time needed to run a simulation. 

In this particular case, since Ry = —A x £g> 1 — 1 <S> A t is the sum of 
two tensorized operators which commute with one another, we can use an 
approach described in [11] . This method consists in using an approximate 
expansion of the inverse of the Laplacian operator, constructed as follows. 
The function h : x G [xq, +oo) i— > - (where xq is a positive real number) can 
be approximated by a sum of exponential functions of the form 



1 N 

x ^ 

i=i 



Qe 



for some JVeN* and where (C;)i</<at and (q)i<ktv a two sets of well-chosen 
real numbers, depending on x$. Provided that xq satisfies xq < min(l, Xf, A'), 
where Af (respectively A') is the lowest eigenvalue of the operator — A x on 
Hq(X) with respect to the L 2 (X) scalar product (respectively the lowest 
eigenvalue of the operator — A t on Hq{1~) with respect to the L 2 (T) scalar 
product), since both the operators — A x <g> 1 and —1 (g) A t commute, Ry 1 can 
be approximated by 



Ry' ~ Eli Qe-^-^® 1 - 1 ®^ 



N^ n - ClAx ^ p - Cl A t ( 26 ) 



The computation of the expansion ( 126|) only involves the computation of 
the exponential of small-dimensional operators. But of course, to have a 
reliable approximation of this operator, the number N of terms in the above 
approximation may be very large. Besides, an explicit expansion is not always 
available for a general operator Ry 1 . 

The algorithms presented in the following sections are attempts to find 
numerical methods which 

• avoid the conditioning problem inherent to the method described in 
Section 13. H 
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• avoid the use of inverse operators such as R v in the approach using 
the dual norm. 

Of course, a natural idea would be to find a suitable norm to minimize 
the residual to avoid the conditioning and inversion problems. So far, no 
norms with such properties have been proposed. 

In Section HI we present algorithms already existing in the literature, 
namely those suggested by Anthony Nouy [TB] and Alexe'i Lozinski [14] . In 
Section \5\ a new algorithm is proposed. The known partial convergence 
results for these methods are presented and the numerical implementation of 
the algorithms are detailed. 

4 Algorithms based on dual formulations 

In this section, we present some classes of algorithms based on dual formu- 
lations of the non-symmetric problem (TT6j) . 

4.1 MiniMax algorithm 

A first algorithm based on a dual formulation of problem (TT6j) is the MiniMax 
algorithm proposed by Nouy |16| . 
The algorithm reads as follows: 

1. let u = and n — 1; 

2. let (r n ,r n ,s n ,s n ) G V£ x V t 2 such that 

(r n , r n , s n , s n ) G arg max min J n {r®s,r®s), (27) 

(r,s)eV x xV t (r,s)&V x xV t 

where for all v, v G V, 

J n (v,v) = -\\v\\v - a(u n -i + v,v) + l(v); 

3. set u n = u n -i + r n <g) s n and n = n + 1. 

At each iteration n G N*, the computation of a quadruplet (r n , r n , s n , !s n ) G 
V% x V t 2 satisfying (127|) is done by solving the stationarity equations 

f a(u n -i + r n ® s n ,r n ® 5s + 5r® s n ) = l(r n ® 5s + 5r ® s n ), V(5f, 5s) G V x x V u 
\ a(r n ®5s + 5r® s n , r n ® s n ) = (r n <g> 5s + 5r <g> s n , r n ® s n ) v , V(<5r, 5s) G V x x V*. 

(28) 
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In practice, for each n6f, these equations are solved through a fixed- 
point procedure where the pairs (r n ,r n ) G V£ and (s n , s n ) G V^ 2 are computed 
iteratively. More precisely, the fixed-point algorithm reads: 

• set m = 0, and choose an inital guess (r n \rn \ Sn\s^j G V x x V t 2 ; 

• find (ri m+1) ,?i m+1) ) G such that 

a ^ n _! + r n m+1) ® 4 m) , 5r ® 3< n m) ) = I (ft® Si m) ) , V<5f G 

a («5r ® S l m) , ?l m+1) ® Sl m) ) = (<5r ® ^ m) ,r^ +1) ® ^ m) )^ , V5r G 

• find (^ m+1) ,3l m+1) ) G if such that 

a (u n _ x + ri m+1) ® , ?< m+1) (g) 5?) = Z (?< m+1) ® <5s) , G K , 

a \r n m+l) ® 5s, 7i m+1) ® = (r n m+l) ® to, ri m+1) ® S ( m+1) )^ , Vto G Vfc 

• set m = m + 1 . 

In [IT], it is proved that in the case when a = a x ® a t where a x is a 
continuous bilinear form on 14 x V,. and a t a continuous bilinear form on 
V t x V t and 1^ = 14®^, the algorithm converges. However, there is no 
convergence result in the full general case. 

4.2 Greedy algorithms for Banach spaces 

Another family of dual greedy algorithms is inspired from the methods sug- 
gested by Temlyakov in [20] for Banach spaces and was proposed by Lozin- 
ski [13] in order to deal with the resolution of high-dimensional problems of 
the form (TIT)]) . 



4.2.1 Greedy algorithms for general Banach spaces 

For the sake of simplicity, let us present two particular greedy algorithms 
proposed by Temlyakov in the context of Banach spaces, namely the X- 
Greedy and the Dual Greedy algorithms. 
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Let (X, || • || x) be a reflexive Banach space and T> a dictionary of X, i.e. 

a subset of X such that for all g G T>, \\g\\x = 1 and Span(P) = X. Let us 
also denote by X* the dual space of X. 

Let / G X. The aim of both the Dual Greedy and the X-Greedy algo- 
rithms is to give an approximation of / as a linear combination of vectors of 
the dictionary T>. These numerical methods are generalizations of the Pure 
Greedy algorithm, which is defined for Hilbert spaces. When X is a Hilbert 
space endowed with the scalar product (-,-)x, the Pure Greedy algorithm 
can be interpreted in two equivalent ways, namely: 

Pure Greedy algorithm (1): 

1. let /o = 0, ro = / and n = 1; 

2. let g n ET> and a n G K such that (assuming existence) 

||?V-i - a n g n \\ x = min ||r n _i - a#||x; 

3. let /„ = / n _i + r n = r n _i - a n g n and n = n+ 1; 
and 

Pure Greedy algorithm (2): 

1. let /o = 0, r = / and n = 1; 

2. let g n ^T> such that (assuming existence) 

(r n -i,9n)x = m&x(r n _ 1 ,g) x ; 
gev 

3. let a n G M such that 

||r n _i - Qf n ^n||x = min ||r n _i - ag n \\ x ; 

4. let / n = / n _i + a„5( n , r n = r n _ x - a n g n and n = n+ 1. 

When X is a Hilbert space, the two versions of the Pure Greedy algorithm 
are equivalent, but this is not the case anymore as soon as X is a general 
Banach space. 

The X-Greedy algorithm corresponds to the extension of the first version 
of the Pure Greedy algorithm: 
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1. let /o = 0, ro = / and n = 1; 

2. let g n aT> and a n G R such that (assuming existence) 

||r n _i - a n ^ n || x = min ||r n _i - a#|| x ; (29) 

3. let /„ = / n _i + a n g n , r n = r n _i - a n g n and n = n+l. 

The .Dim/ Greedy algorithm generalizes the second version of the Pure 
Greedy algorithm and is slightly more subtle. It is based on the notion of 
peak functional. For any non-zero element / G X, we say that Ff G X' is a 
peak functional for / if H-f/Hx* = 1 an d Ff(f) = \\f\\x- The Dual Greedy 
algorithm reads: 

1. let /o = 0, r = / and n — 1; 

2. let F rn _ 1 G X* be a peak functional for r n _! and let g n G £> such that 
(assuming existence) 

g n G argmax F Tn _ t (g) ; (30) 



3. let a n G K such that 

a n G argmin ||r n _i - a#n||x; (31) 



4. let /„ = + r n = r n _! - a n # n and n = n+ 1. 

Slightly modified versions (relaxed versions) of the X-Greedy and Dual 
Greedy algorithms are proved to converge in (20] provided that the space X 
and the dictionary T> satisfy some additional assumptions, detailed below. 

We define the modulus of smoothness of the Banach space X by 
V/? G E, p(P) := sup (h\\x + Py\\ x + \\x - py\\ x ) - 



IMIx=|l3/llx =1 



2 



The Banach space (X, \\ ■ \\x) is said to be uniformly smooth j20] if 

lim^ = 0. 
23 



Let us point out that if a Banach space (X, \\ • \\x) is uniformly smooth, then 
the mapping G:i6l4 ll^llx is Frechet-differentiable. 

The relaxed versions of the X-Greedy and Dual Greedy algorithms are 
proved to converge |2U] provided that 

(Bl) Span(£>)"'" x = X; 

(B2) WD is weakly closed in X; 

(B3) X is a uniformly smooth Banach space. 

We do not write here these relaxed versions of the algorithms for the sake of 
brevity and refer to |2U| . 



4.2.2 Special Banach spaces for non-symmetric high-dimensional 
problems 

Let us now present how these ideas were adapted by Lozinski to the case of 
high-dimensional non-symmetric problems. We begin here with the descrip- 
tion of the particular Banach spaces involved. Let us assume in the rest of 
Section 14721 that the operator A -1 : V — >■ V is bounded. 



A Banach space with good theoretical properties but which can- 
not be used in practice 

The space V is now endowed with the following dual norm 

Vv G V, ||f |U = sup 77 1 - = \\Av\\v = \\Av\\y- 
w£V,w^0 \\w\\v 

Actually, since the linear operator A is bounded on V, the space (V, || ■ 
is a reflexive Banach space whose dual space is (V, \\ ■ ||(^.*)-i) where 

V« e V, H^.j-i = sup = WiAT^Wv- 

wev,w^o \\ w \\v 

Let us show that the Banach space (V, || • m) and the dictionary 
V = {r ® s, r € V x , s G V t , \\r ® s\\ A = 1} 
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satisfy assumptions (Bl), (B2) and (B3). 

Let us begin with the proof of (Bl) and (B2). Since the set of tensor 
product functions 



£ = { r <g> s, r G V x , s E V t } = WD 



-(V.IMIv-) 



is assumed to be weakly closed in (V, || ■ and to satisfy Span(S) ' = V 
(assumptions (Al) and (A2)), (Bl) and (B2) are direct consequences of the 
fact that A and A~ l belong to the space £(V, V) (i.e. are bounded operators). 
For instance, £((V, || • ||y),R) = £ ((V, || • |U),R) since 

1 

V7 G £((V,|| • ||y),M) , tt-t z |KIU((v,|H| v ),iR) < IKIU((v,||-|U),k) < ||^- _1 |U(i/,y)IKIl£((y,||-||v),: 

||^||£(V,V) 

Let us now prove (B3). Since the operator A is invertible, the modulus 
p_4 of smoothness of (V, \\ ■ is equal to the modulus of smoothness p of 
(V, || • || v ). Indeed, for all /3 G M, 

Pa(P)= sup (h\\v + Pw\\ A +\\v-Pw\\A)-l 

= sup (-(\\Av + 0Aw\\ v +\\Av-PAw\\ v )-l 

V,W£V, \\Av\\y = \\Aw\\y = l \^ 

= sup ( -(\\v + )3w\\v + \\v - /3w\\v) - 1 

v,w£V, \\v || v" = ||tu|| y=l V 
= Ptf). 



Since (V, || • ||y) is a Hilbert space, 

and (V", || ■ m) is a uniformly smooth Banach space. 

To implement the X-Greedy or Dual Greedy algorithms in practice in this 
context, one needs to compute the norm || • \\ A (see (129]) and (l3Tj) ). Since for 
all v E V, \\v\\a = \\Av\\v = WR^AvWv, this requires the resolution of several 
intermediate low- or high-dimensional problems to compute the inverse of the 
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operator Ry. The same issues as those described in Section 1X21 have to be 
faced. 



Another more practical Banach space 

The idea of Lozinski is to replace this norm by a weaker one, easier to 
compute, 

w r- t/ II II a(v,r®s) 

Vv G V, \\v\\ iA = sup — . (32) 

(r,s)€V x xV t , r®s^0 \\ r & S \\V 

Actually, denoting by || • ||j the injective norm on V j8], defined by 
w c t/ ll II (v,r®s) v 

VV G V, \\v\\i = SUp — ■ r: , 

(r,s)eV x XV t , r®s^0 || r ® s||y 

it holds that for all ueV, IMIc4 — ll^^lli- Reasoning as above, the Banach 
space (V, || • || 1,4) has exactly the same properties as (V, || ■ ||j). 

Since S is weakly closed in V, and since for all v G V, ||v||j < |H|v, S 
is also weakly closed in (V, || • ||j). But, in the full general case, the Banach 
space (V, || ■ || i) and hence the Banach space (V, || • ||^) are not uniformly 
smooth. Actually, these spaces may not even be reflexive. Indeed, let us 
assume that V = V x <g> V t and that || ■ ||y is the associated cross-norm, in 
other words that for all (r,s) G V x x Vt, ||r <E> s\\v = \\ r \\v x \\ s \\vf ^ holds 
that [S] (V, || • ||i) is isomorphous to K,(V x ,V t ), the Banach space of compact 
operators from V x to Vt endowed with the operator norm. Since JC(V x ,Vt) 
is not a reflexive space (JC(V x ,V t )* = 61^, V*) and &i(V x , V t )* = 2{V x ,V t ) 
where &i(V x , Vt) denotes the set of trace-class operators from V x to Vt), there 
is no guarantee of convergence of the relaxed versions of the X-Greedy or Dual 
Greedy algorithms presented above. 



The finite-dimensional cross-norm case 

However, in the case when V x and V t are finite-dimensional and V = 
V x ®Vt, the spaces K.(V X , V t ) and £(V X , V t ) are identical. The space (V, || • ||j) 
is then reflexive and uniformly smooth. Indeed, if V x = M. mx and V t = M mt , 
since || • \\y x and || • ||y t both derive from the scalar products (•, -)y x and (•, -)y t , 
there exist two invertible matrices P E M. mxXmx and Q £ M miXmi such that 

Vr G M. mx , \\r\\ Vx = \\Pr\\ Fmx , 
VsGM mi , \\s\\ Vt = \\Qs\\ Fm . 



2(3 



where || • \\F mx and || ■ \\F mt denote respectively the Frobenius norms on K™ 1 
and M m *. Thus, (V, \\ ■ ||j) is isometrically isomorphic to W nxXmt seen as 
fc(V x , Vt), endowed with the norm 

\\QMr\\ F 

VM G M m - Xm *, \\M\\i = sup 7 " mt = \\QMP-%, 

where 

\\Mr\\ F 

MM G ]R m * xm ', ||M|| 2 = sup 



reR m x,rj!=0 \\ r \\F, 



■m x 



Actually, (IR ma:Xmt , || • || 2 ) is a uniformly smooth Banach space [19J. Thus, 
greedy algorithms for Banach spaces do converge in this setting. 



4.2.3 Practical implementation of the algorithms 

X-Greedy algorithm 

The X-Greedy algorithm reads: 

1. let uq = and n — 1; 

2. find (r n , s n ) G V x x V t such that 

{r n ,s n )£ argmin \\u - u n -i - r <g) s\\u] (33) 

(r,s)&V x xV t 

3. set u n = « n _i + r n ® s n and n — n + 1. 

From the definition of the norm || ■ ||j_4 (see (132]) ). the second step of the 
algorithm can be rewritten as 

(2) find (r n , s n ) G V x x Vt such that 

(r n ,s n j G argmm (rs)el4x y t sup (? r^ e y ;cX y t pgfj^ 

= argmin (7 , s)eV:cXyt sup^^^ ¥ ®S\\ y =\ l(r®a)- aK-i +r®s,r®s). 

(34) 
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From a practical point of view, at each iteration n G W, the functions 
{ r n,r n , s n ,H n ) G V 2 x V t 2 are obtained by solving the stationarity equations 
associated with ( 134")) . namely by solving the following coupled problem 

( find (r n ,r n , s n , s n ) G V 2 x V t 2 such that for all (5r,5r,5s,5J) G V 2 x V^ 2 , 

Cg> Ss + 5r ® s n )y + a{u n -i + ?"n ® s n , r n <E> 5s + 5r <g> s n ) — l(r n <g> 5s + 5r <E> s n )i 
[ a(r„ (g) 5s + <5r <g> s n , f n <S> s n ) = 0. 

The X-Greedy algorithm has not been implemented in practice yet. In- 
deed, it is not clear how to compute a solution of the above stationarity 
equations since using a fixed-point algorithm procedure similar to the one 
presented in Section 14.11 for the MiniMax algorithm would always lead to 
T n ® s n = 0, due to the form of the second equation. 

Dual Greedy algorithm 

Let us describe here how Lozinski adapted the Dual Greedy algorithm for 
the resolution of high-dimensional non-symmetric linear problems. 

A remaining issue concerns the construction of a peak functional F Tn _ 1 
for the residual r n _i = u — w n _i which is used in the second step of the 
algorithm (I5U|) . Actually, the true peak functional is not computed but only 
an approximation of this functional by an optimal tensor product function 
in a sense which is made precise below. 

The adapted Dual Greedy algorithm reads: 

1. set Mo = and n = 1; 

2. (computation of an approximate peak functional for the residual r„_i = 
u — u n -\ with a tensor product function) find (r n ,'s n ) G V x x Vt such 
that 

(^n,s n ) G argmax a{u—u n -\,r®s)= argmax l{r®s)—a{u n -\,r®s)] 

(r s s)eV x xVt, ||rig>s1|y=l {r,s)eV x xVt s ||f®Sl|y=l 

(35) 

3. (Step 2 of the Dual Greedy algorithm, see (I3"U|) ) find (r n , s n ) G V x x V t 
such that 

(r n , s n ) G argmax a(r <g> s, r„ ® s n ); (36) 

(r,s)£V x xV t , 1^0511^=1 
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4. (Step 3 of the Dual Greedy algorithm, see (131 j) ) find a n G 1R such that 

a n G argmin \\u—u n -i—ar n &)s n \\i_A = argmin sup /(r®5) — a{u n _i+ar n ®s n) r®s 

(37) 

5. set u n = u n -i + a n r n ® s n and n = n + 1. 

The Euler-Lagrange equations associated with the minimization problem 
d35D read: for all (8r, 5s)EV x xV t , 

K(r n ®Sn,r n ®5s + 5r(3s n } v = l(r n (g)5s + 5r(g)s n )-a(u n _ 1 ,r n ®5s + 5r®s n ), 

for some A n G 1R satisfying AnH^fgJ^Hy = A„ = l(r n <^)'s n ) — a(ti n ,_ 1; r n ® 7 n ). 

The Euler-Lagrange equations associated to the minimization problem 
( |36|) can be rewritten as follows: (r n , s n ) G 14. x V t is solution of 

a(r n ®5s + 5r®s n ,r n ®s n ) = n n a(r n ®5s + 5r®s n ,r n ®'s n ), V(8r,8s) G V x xV u 

where (r„,s"„) G 14 x is such that 

(f„,s n ) G argmax a(r ra <g) s n , r <g) s), (38) 

(r,s)6l4xVt, ||fi»a||v=l 

and /i n = a(r n ® s n ,r n ® s„). Besides, if the pair (r^,^,) satisfies ( 13"B"j) . it 
holds that 

a(r n (gis n ,f n ®<5s+<Sf<gis n ) = fcVi(^n®^n>^»®^"+<^®sii)v> V(Sr,$s) G 14x14, 

with z/„ = a(r n (g) s n ,r n (g) s~ n ) = ||r n (8) s^U^ = 1. This yields to a coupled 
problem on (r n ,s n ) and (r n ,s~ n ). Lozinski the noticed that, if (r n ,s n ) is 
solution to ( )36|) . (f^s^) = (r n ,s n ) is solution of ( 138]) in the sense that 

(f n , s n ) G argmax a(r n <g) s n , r <g) s). 

(r,s)&V x xV t , \\r(gis\\ v =l 

The Euler-Lagrange equations can be rewritten as 

a(r n ®s n , r n ®6s + 5r®s n ) = (r n ®s n ,r n (g)5s + 5r®s n ) v , \/(Sr,8s) G V x xV t , 
by taking /i n = a(r n <g> s n ,r n <g> s n ) = \\r n <g> s n \\ iA = 1. 
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The Euler equations associated with (I37p read 

a(u - w n _i - a n r n ® s n , r n <g> s n ) = 0, 

yielding to a n = a{u — u n -i,r n (g) s" n ) = A„. 

Finally, replacing a n r n ® s ra by r n ® s„, the iterations of the Dual Greedy 
algorithm are computed in practice as follows: 

1. let uq = and n = 1; 

2. find (r n , r n , s n ,s n ) G V^ 2 x VJ 2 such that for all (5r, 5r, 5s, 57) G V£ x V^ 2 , 

f (r n <g> s n ,r n ® 5s + 5r <g> 7 n ) v = l(r n <8> <5s + 5r <g> s n ) - a(« n _i,r n (g> 5s + <5f <g> s n ), 
\ a(r n <g) s n , r n <g> 57 + 5r®7 n ) = l(r n <g> 57 + 5r ®7 n ) - a(u n _ 1 ,r n <g) 57 + 5r ® s n ); 

3. set u n = w n _i + r n <g> s n and n — n + 1. 

These equations lead to two decoupled problems on (r n , s n ) and (r n , 7 n ). Each 
of them is solved through a fixed-point procedure similar to (TT5]) described 
in details in Section 16.3.11 Some numerical tests are presented in Section [61 
which illustrate the convergence of this algorithm. 

5 The Decomposition algorithm 

Let us now present a new algorithm based on a decomposition of the bilinear 
form a. The bilinear form a can always be written as 

a(;-) = b a (; •)+&(•,•) (39) 

where b s (-, •) is a symmetric coercive continuous bilinear form on V x 1/ and 
&(•, •) a (not necessarily symmetric) continuous bilinear form on V x V. In 
the sequel, b s (-, •) (respectively &(•,•)) will be refered to as the implicit (re- 
spectively explicit) part of a(-, ■). This terminology will be explained below. 

Example 5.1. Let us introduce a s (respectively a as ) the symmetric part (re- 
spectively antisymmetric part) ofa(-,-), defined by: 

Wv, w G V, a s (v, w ) — — (a(v, w) + a(w, v)) , (40) 
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and 



\fv, w G V, a as (v, w) = - (a(v, w) - a(w, v)) . 



(41) 



Provided that a s is coercive, the decomposition a(-,-) = a s (-,-) + a as (-,-) is 
admissible in the sense of f3W . 

Let us denote by B s and B the bounded operators on V defined by 



Since b s is coercive, the operator B s is invertible. 

The principle of the algorithm we consider in this section to solve problem 
(TT6]) consists in expliciting the part b of the bilinear form as a right-hand 
side source term. More precisely, one can consider the following fixed-point 
algorithm: 

1. choose a starting guess uq G V and set n = 1; 

2. let u n be the unique solution to 



3. set n — n + 1. 

In other words, u n = F(it re _i) where for all u 6 V, F(v) = B~ l (C — Bv). 
If ||i3,7 1 i3||£(yy) < 1, the mapping F : V — >• V is a contraction and it fol- 
lows from the Picard fixed-point theorem that the sequence (u n ) n6 N strongly 
converges in V towards the solution u of the initial problem (I16jl . 

A natural approach thus consists in solving problem (142)) at each itera- 
tion n G N* using a standard greedy procedure. Provided that the greedy 
expansion obtained at each iteration n G N* is accurate enough, the sequence 
{ u n)neN given by this algorithm strongly converges in V towards u. 

However, the principle of this method requires to compute a full greedy 
loop at each iteration n G N*. In order to save computational time, we now 
introduce the following algorithm, in which only one tensor product function 
is computed at each iteration n G N*. 

Decomposition algorithm: 



Wv,w G V, b s (v,w) = (B s v,w) v , 
Wv,w G V, b(v,w) = (Bv,w)v- 



Wv G V, b s (u n ,v) = l(v) - b(u n -i,v); 



(42) 



31 



1. let % = and n — 1; 

2. find (r n , s n ) E V x x V t such that 

(r n ,s n )e argmin ~b s {u n -x+r®s,u n - 1 +r®s)-l(r®s)-b(u n - 1 ,r®s)] 

{r,s)eV x xV t 2 

(43) 

3. set u n = u n -i + r n (g) s n and set n = n + 1. 

The bilinear form &(•, •) is explicited as a right-hand side, whereas 
remains implicit. This justifies the terminology introduced in the beginning 
of the section. 

Equation ( 143]) can be rewritten equivalently as 

(r n , s n ) G argmin (r tS ) eVxXVt \b s (r <g> s,r ® s) - l(r ® s) - 6 s (w„_i + r ® s) - &(w„_i, r <g> s) 
= argmin (rs)6V3;XVt \b s (r ® s,r <8> s) - Z(r g> s) - a(u n _i, r <g> s). 

This means that for each n G N*, tensor product solution to 

the first iteration of the greedy algorithm applied to the following symmetric 
coercive problem 

{find u G V such that 
W G V, b s (w, v) = l(v) - i> s (u n -i, u) - &(«n-i, = - a(w n _i, u). 

As explained above, such an algorithm is expected to converge only if the 
norm ||Sj 1 jB||£(y ) y) is small enough. Actually, in the case when the spaces 
V x and V t are finite dimensional, the following result holds: 

Proposition 5.1. If V x and Vt are finite dimensional, there exists k > 
such that if 

\\B^B\\ z{ y y) < k, (44) 

then the sequence (u n ) ne ^* defined by the Decomposition algorithm strongly 
converges to u in V . 

Proof. Since V t and V x are assumed to be finite dimensional, from assumption 
(Al), so is V. Let k := HBj^l^^y)- Let us denote by (•, -)y = b s (-, •) the 
scalar product on V induced by the symmetric bilinear form b s , and by || • ||y 
the associated norm. 
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Let £ G V and B G £(V, 7) be denned by 

Wv e V, l(v) = (£^v)y, 
Wv,w G V, b(v,w) = (Bv,w)y. 



Actually, B = B S 1 B and C = B s 1 £. This implies that for all v G V, 
\\Av\\ v < n\\v\\y. Besides, I + B = B~ 1 (B S + B) = B~ 1 A, where X denotes 
the identity operator on V. 

The proof relies on the fact that all norms are equivalent in finite dimen- 
sion. In particular, let us define the injective norm || • ||rby 



Wv G V, \\v\\z= sup 



(v, r <g> s) v 



(r, s )ev x xv t \\ r ®s\\v 
Then, there exists a > such that 

Wv G V, \\v\\ v < ct\\v\\z. 

Besides, for all v G V, |H|r< \\ v \\v- 

Let us introduce U n := C — Bu n — u n = B~ l (C — Au n ). For all n G N*, 
U n is the vector of V such that 

Wv G V, (U n ,v) v = l(v) -a(u n ,v). 

For n > 1, the Euler equations associated to ( |43|) read: 

& s (ti„_i+r„®s n ,r n ®<fe+5r<g)s n ) = l(r n ®5s+5r<gis n )—b(u n -i,r n <8>8s+5r<gis n ), W(5r,5s) G V x xV t . 

As a consequence, l(r n (g) s n ) - a(u n ^i z r n <g> s n ) -b s (r n <g> s n ) = — r n (8) 

s n ,r n ®s n ) v = 0, which implies that ||Z4_i||~ = ||W n -i-r n <g)s ri ||~+||r ri <8)s w ||| f . 
Furthermore, from Lemma we have 

ii ii (U n -i,r®s)y ~ 

F„, (g) s n II v = sup — = ||64_i||r- 

(r,s)eVixVt, r®s/0 ll r ® s lly 
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It holds that 



fe-i||~ - \\U n \\ 2 y = \\&n-i\\ 2 y ~ \\Un-i ~ r n <g> s n - Br n ® s n \\ 2 ~ 

= - \\U n -i - r n ® s n \\ v - \\Br n ® s n ||| 

+2(Br n <g> s n ,U n -i - r n ® s n ) v 

= \\r n <8> s n \\~. - \\Br n <g> s n \\~ - 2(Br n <g> s n ,W n _i - r n <g> s n ) v 

> (1 - n 2 )\\r n <g> s n ||~ - 2«||r n <g> s n || v ||Z4_i - r n <g> s n || v 

> (1 - K 2 )\\r n <g> s n \\ v - 2n\\r n ® s n \\y\\U n -i\\y 

> (1 - /t 2 )||r n O s n ||~ - 2«:a||r n <g> sj v ||Z4_i||r 
= (1 - n 2 - 2a.K)\\r n <g> s n ||~ 

= (1-k 2 -2aK)||W n _i |||. 



If « is small enough to ensure that 1 — n 2 — 2an > 0, the sequence (||W n ||^) ng N* 
is non-increasing, hence convergent. Thus, the series of general terms ( ||£/ n _i || ~— 
l|W n ||~)neN* and (||W n ||~) neN . are convergent. This yields 

|ft||r— ^0, 

n— >oo 

and, as all norms are equivalent in finite dimension, 

felly — >«■ 

n— >oo 

If K < 1, the operator (X + S)" 1 is continuous from V to V for the norm 
|| ■ ||y, and 

||n-^Hv = \\A- l C-u n \\ v = ^(B^A)' 1 B^C - Un|_ = ||(X+^)- 1 (£-(X+^) Mn 

Since the norm || • || y is equivalent to the norm || • ||y, we obtain the desired 
result. □ 

Unfortunately, the rate k in (144p strongly depends on the dimensions of 
14 and Vj, as shown in the proof of Proposition 15. II However, in some simple 
numerical experiments we performed with this algorithm, the rate k does not 
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seem to depend on the dimension, as illustrated in Section The analysis 
of this numerical observation is work in progress. 

A possible way to obtain a decomposition such as ( 131?|) satisfying f l4"4")) is 
to consider preconditioners for the initial antisymmetric problem. Problem 
( TLB"]) can be rewritten as 



where C is a well-chosen continuous linear operator on V such that ||X — 
C»4||£(v,y) is as small as possible. With this formulation, one can consider 
the following algorithm: 

1. let u = and n — 1; 

2. find (r n , s n ) G V x x V t such that 

(r n , s n ) e argmin ^ (u n _i + r <g> s, u n _i + r <gi s) y -((X-C^l)u„_ 1 +C£, r 



3. set n n = n n _i + r n (g) s„ and n — n + 1. 

The bilinear form associated with the continuous operator C^4 needs to be 
easy to compute on tensor product functions. At least in the case when V is 
finite dimensional, if ||X — C^4.|| £(v,v) is small enough, the sequence (ti n )„ e N* 
converges strongly in V towards the solution u of (fT6l) . However, finding a 
suitable operator C is not an easy task in general. 

6 Numerical results 

6.1 Presentation of the toy problems 

In this section, we present the two toy problems we consider in these numer- 
ical tests. Let b x , b t G M, m x = m t = 1 and X = T = (—1, 1). 

The first toy problem is inspired from Example 12.11 In the continuous 
setting, we define 



. V:=LL V ((-l,l),m cr ((-l,l),C)) = L2 ((-l,l),C)®m cr ((-l,l),C); 




(45) 



(46) 
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14 := ^ er ((-l,l),C); 
^:=L^((-1,1),C); 

/e^ er ((-l,l) 2 ,C); 

for all u, v G V, 



•i r i 

a(u, 



J I (V x u ■ V x v + (b x ■ V x u)v + uv) , 



and 

»i pi 



l{v) := J J Jv. 

The second toy problem is the following 

• ^:=^ er ((-l ; l) 2 ,C)); 
. V x :=Hl ex {{-lMC); 
. y,:=^ cr ((-l,l),C); 

. feLl ei ((-l,iy,C); 

• for all u, v G V, 

a(u, v) := J J (V x u ■ V x v + V t u ■ V t v) 

/ ((b x ■ V x u + b t ■ V t u)v + uv) 
1 J-i 



and 

l(v):= fjjv. 

It can be easily checked in the two cases that the Hilbert spaces V, V x 
and Vt satisfy assumptions (Al) and (A2). Besides, the sesquilinear form 
a : V x V — > C is continuous and satisfies assumption (A3). Our aim is to 
approximate the function u G V such that 

G V, a(u,v) = l(v). (47) 
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We introduce finite-dimensional discretization spaces defined as follows. 
For all k G Z, let e x k : X 3 x ^ -±e ikx and e\ : T 3 t ^ 4^e ikt . For 
N x ,N t eW, we define 

V?*:=Spaa{e%, -N x < k < N x } , 

^*:=Span{4 —N t <k< N t } , 

and V^-'^ := Vf» ® Vf J . 

We then consider the Galerkin approximation of problem f|47j) in the 
discretization space V Nx ' Nt , i.e. find u G V Nx ' Nt such that 

Vt; G V^"-^, a(u,v) = l{y). (48) 

For u solution of (HS). we denote by U = (u kl ) ]kl < Nx , mNt G C( 2JV * +1 ) x ( 2JVt+1 ) 
the matrix such that 

N x N t 

u= Yl E u klzl®Zu 

k=-N x l=-N t 

and by F = (f k i) ]k \< Nx , \i\<N t e & 2N * +1 ^ 2N ^ the matrix such that 

N x N t 

f= E E /^® e i- 

In this discrete setting, the first problem is equivalent to: find U G £,( 2N ^+ 1 ) x ( 2N t+ 1 ) 
such that 

D x U + b x N x U + U = F 
and the second problem equivalent to: find U G c( 2JV *+ 1 ) x ( 2JV t+ 1 ) such that 

D X U + b x N x U + UD t + btUN t + U = F (49) 

where D x , N x G C( 2JV * +1 ) X ( 2JV * +1 ), D t , N t G 2Nt+1 ^ 2Nt+l \ and for all -N x < 
k, k' < N x and —N t <l,l' < N t , 

{D x )kk' = \k\ 2 Skk', 

{D t )u> = |/| 2 ^', 

(N x ) kk > = ik5 kk >, 

{N t ) lv = ilSw. 
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6.2 Tests with the Decomposition algorithm 

Let us begin with the numerical tests performed to illustrate the convergence 
of the Decomposition algorithm presented in Section [5j 

6.2.1 The fixed-point loop 

Before presenting the numerical results obtained with this algorithm, we 
would like to discuss the fixed-point procedure used in practice in order to 
compute the pair of functions (r n , s n ) G V t Nt x V^ x , solution of (143)) at each 
iteration n G N* of the Decomposition algorithm. 
The algorithm reads as follows: 

• choose (Vi 0) , s ( n ] ) G V t Nt x V*' and set m = 1; 

• find (r [ ™\st ] ) G ^ x V?' such that 

6 S (rt ] ® 4 m_1) , 5r ® st~ 1] ) = I (Sr ® - a ( Un _ 1? 5r ® s^) , V<5r 

6 S (ri m) ® ,i m) , ri m) s Ss) = I (ri m) ® Ss) - a (u^-i, ri m) ® 5s) , V5s G V t m ; 

• set m = m + 1. 

This fixed-point procedure exhibits of the exponential convergence rate 
which is numerically observed for standard greedy algorithms for symmetric 
coercive problems, as discussed in Section [TJ The stopping criterion we choose 
for simulations presented in Section 16.2.21 is the following: \\rn £g> — 
?"i m-1 ' ) ® si m_1 ^|| < e where || • || denotes the Frobenius norm and e = 10 -8 . 

6.2.2 Numerical results 

We consider the second problem presented in Section [6TT| where / G Lp er ((— 1, l) 2 , C) 
is chosen such that 

kez zez 

with f H = |fc|2+ ^ |2+1 for all k, I G Z. 

We decompose the bilinear form a(-, •) as a(-, ■) = b s (-, ■) + b(-, ■) where 
b s (-, •) = a s (-, •) is the symmetric part of a(-, •) defined in ( l40l) and b(-, •) = 
a as (-, •) is the antisymmetric part of a(-, •) defined in (14ip . 
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Let us recall that Proposition 15.11 states that, in the finite-dimensional 
case, there exists a rate n (which depends on the dimension of the Hilbert 
spaces) small enough such that if H^" 1 ^ £(v,v) < the Decomposition 
algorithm is ensured to converge. In the numerical simulations presented 
below, we have witnessed that there exists a threshold rate k such that 

• if ||i37 1 J B||£(\/,y) < k, the algorithm converges; 

• if = the algorithm does not converge, but the norm of 
the residual remains bounded; 

• if > k, the algorithm does not converge, and the norm of 
the residual blows up. 

Besides, the rate k seems to be independent of the dimension of the Hilbert 
spaces. 

For nEW, let U n G C {2Nx+1)x{2Nt+1) denote the approximation of U 
solution of ( |49|) given by the algorithm at the n th iteration. The following 
three figures show the evolution of the logarithm of the norm of the residual 
|| F - D x U n - U n D t - b x N x U n - h t U n N t - U\\ &2 (where || • ||e 2 denotes the 
Frobenius norm), for different values of N = N x = N t and b = b x = b t as a 
function of n. We observe numerically that in this case, the limiting rate is 
obtained for b = 1.5 for any value of N. 

We also performed another set of tests with a modified version of the 
Decomposition algorithm, in order to increase the threshold rate k, in a sense 
which will be precised below. The modified Decomposition algorithm reads 
as follows for a G [1, +oo): 

1. let uq = and n = 1; 

2. find (r n , s n ) G V x x V t such that 

Oi 

( r mS n )G argmin — b s (r (8) s, r (g> s) — l(r ® s) — a(it n _i, r (g> s); (50) 

(r,s)eVxXV t 2 

3. set u n = u n ^i + r n ® s n and set n — n + 1. 

Let us point out that this algorithms is equivalent to the standard Decom- 
position algorithm presented in Section in the case when a — 1. 
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Figure 2: Evolution of log 10 (||F - D x U n - U n D t - b x N x U n - b t U n N t - U\\ 6a ) 
as a function of n for = 20 and different values of b. 
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Figure 3: Evolution of log 10 (||F - D x U n - U n D t - b x N x U n - b t U n N t - U\\ 6a ) 
as a function of n for = 50 and different values of b. 
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Equivalently, for each n G N*, r n ® s n is a tensor product solution to 
the first iteration of the greedy algorithm applied to the symmetric coercive 
problem 

J find u G V such that 

|Vi> G V, ab s (u,v) = l(v) — a(u n -i,v). 

We observe that this algorithm has the same convergence properties as 
the standard Decomposition algorithm, i.e. there exists a threshold rate 
n a = WB^BW^v) such that 

• if ||i57 1 i3||£(y > y) < k q , the algorithm converges; 

• if 1 1 i?^ ^ 1 1 mv,v) — K a , the algorithm does not converge, but the norm 
of the residual remains bounded; 

• if ||i57 1 ,B||£(yy) > n a , the algorithm does not converge, and the norm 
of the residual blows up. 

The rate K a also seems not to depend on the dimension of the Hilbert space. 
Besides, a G [1, +00) 1— > K a seems to be an increasing function. Thus, choos- 
ing a larger value of a seems to lead to an algorithm which is convergent for 
larger values of b. However, the larger a, the smaller the rate of convergence 
of the algorithm for a given value of TV and b. 

The figure below presents the evolution of the logarithm of the norm of 
the residual \\F — D x U n — U n D t — b x N x U n — b t U n N t — U\\q 2 for the second 
problem for a = 2, N = 50 and different values of b. The threshold value of 
b seems to be in this case b = 2.6. 

6.3 Comparison of Galerkin, Minimax, Dual Greedy 
and Decomposition algorithms 

In this section, we present various numerical tests performed to compare the 
performances of the Galerkin, Minimax, Dual Greedy and Decomposition 
algorithms. 

6.3.1 Fixed point procedures 

Before presenting the simulations done to compare the performance of the 
four algorithms, we will detail more precisely the fixed-point procedures used 
for the Galerkin, Dual and MiniMax algorithms. 
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Figure 5: Evolution of log 10 (\\F - D x U n - U n D t - b x N x U n - b t U n N t - U\\ &2 ) 
as a function of n for the modified Decomposition algorithm with a — 2, 
N = 50 and different values of b. 
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Let us first present the fixed-point procedure used in the Galerkin algo- 
rithm. For all n G N*, the method used to compute a pair (r n , s n ) G V x x V t 
solution to ( !T9|) is the following: 



choose (rn \ Sn" 1 ) G V x x V t and set m — 1; 
find ^ri m) , s« m ^ G 14 x F t such that 

a (Vi m) ® 5r ® s^ 1] ) =l(8r® st~ 1] ) - a (u n _ u 6r ® , V5r G ^ 

a (ri m) s s ( n m \ ri m) ® fo) = Z (r< m) ® to) - a (u n ^, ® 5s) , yds G V t Nt ; 



• set m — m + 1. 

All the iterations of the fixed-point procedure are well-defined. However, 
we have seen from Example 12.11 that there are cases where any solution 

(r n ® s n ) G 14 x y t of 

V(<5r, <5s) G V x xV t , a(r n ®s n ,r n (g)5s+5r(g)s n ) = l(r n ®5s+5r®s n )—a(u n _i,r n <g)5s+5r®s n ), 

satisfies r n ® s n = 0. From Example 12.11 this is the case in particular for 
the first problem presented in Section [6.11 with / G Lp Cr ((— 1, l) 2 , C) being 
chosen such that for almost all x,t G (—1,1), f(x,t) = (j)(x — t) with <p G 
Lp er ((— 1, 1), C) a real- valued odd function. The function / = J2k i& fki e %® e \ 
with fki = Si^S-ij + <5_i,fc5i,i for all k, I G Z is an example of such a function. 
On this particular example, we observe numerically that, for n — 1, the 
fixed-point procedure presented above does not converge. 

The practical implementation of the Galerkin algorithm requires the use 
of another stopping criterion as the one used in Section 16.2.11 In the numer- 
ical simulations presented in Section 16. 3. 2} and in the fixed-point procedures 
used in all the other algorithms implemented, the stopping criterion will be 
m < m max with = 20. 

The fixed-point algorithm used for the Decomposition algorithm has been 
detailed in Section I6.2.1[ and the one used for the MiniMax algorithm has 
been described in Section I4TT1 The fixed-point procedure for the Dual Greedy 
algorithm reads as follows: 

• choose ( rii\ si°' ) ) G V x x Vt and set m = 1; 
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while m < m max , find (rh^, si" 1 ' 1 j G V x x V t such that for all (<5r, 6s) G 



x V t , 

^ n m) ®st l - l \57®lt- 1) )v = liSr^lt-^-aiun^S?®^), 
(rl m) g) si m) , rl m) <g <5s)y = l{rk m) (g <5s) - a(w n _i,ri m) <g 5s); 

set m = m + 1 ; 

if 777 > n^maxi Set = and S n = Sn ^ 

choose (^"n \ E V x x Vt and set m — 1; 

while m < m max , find ^ri m ' ) , s n m ' ) J E V x xV t such that for all (<5r, 6s) G 
K x V t , 

a(r„ m) <g Sn , 6r <g s n ) = £((5f <g s n ) - a(if n _i, ^® s„), 
a(ri m) g) s„ , f„ g) (5s) = i(f„ (g 5s) — a(w n _i, f„ <g 5s); 

set m = m + 1 ; 

if m > m max , set r n = r„ m) and s„ = s^. 



6.3.2 Numerical results 

We present here some numerical results obtained for the second problem 
introduced in Section 16.11 with / chosen as in Section 16.2.21 Here iV = 
Nt = N x = 50 and the different figures show the rates of convergence of the 
different algorithms for several values of b = b x = bt- The different values 
of b are the following: 0.01, 0.1, 1 and 2. When b = 2, the Decomposition 
algorithm does not converge. 

We observe numerically that even if the Galerkin algorithm is not well- 
defined through the use of the equations f)19p . using a fixed-point procedure 
as described above and stopping the procedure after a number m max of it- 
erations is enough to observe good convergence properties of the algorithm. 
Besides, the Dual Greedy and MiniMax algorithm are almost as efficient as 
the Galerkin algorithm. 

The Decomposition algorithm is very efficient when the antisymmetric 
part of the bilinear form a(-, ■) is small, but performs badly when b becomes 
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Figure 6: Evolution of log 10 (\\F — D x U n — b x N x U n — U\\@ 2 ) for the different 
algorithms as a function of n with N = 50 and b = 0.01. 
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Figure 7: Evolution of log 10 (\\F — D x U n — b x N x U n — C/||© 2 ) for the different 
algorithms as a function of n with N = 50 and b = 0.1. 
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Figure 8: Evolution of log 10 (\\F — D x U n — b x N x U n — C/||© 2 ) for the different 
algorithms as a function of n with iV = 50 and 6=1. 
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Figure 9: Evolution of log 10 (\\F — D x U n — b x N x U n — C/||© 2 ) for the different 
algorithms as a function of n with N = 50 and 6 = 2. 
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larger. Let us point out that the CPU time needed to compute one tensor 
product in the MiniMax or the Dual Greedy algorithms is twice the time 
needed in the Decomposition or the Galerkin algorithms. Thus, the Decom- 
position algorithm is more efficient than the Dual or the MiniMax algorithm 
for small antisymmetric parts of the bilinear form a(-, •) but, when b becomes 
too large, the algorithm behaves poorly. 

7 Appendix: other algorithms 

In this section, we present two other possible tracks towards the design of 
efficient greedy algorithms for high-dimensional non-symmetric problems, for 
which there is still work in progress. 

7.1 An ill-defined (but converging) algorithm 

In this section, we first present an algorithm whose iterations are ill-defined 
in general but for which we can prove the convergence in the full general 
case. Of course, this algorithm will not be useful in practice but we believe 
the proof is instructive in our context. 

Let a > 0. The algorithm reads as follows: 

1. set uq = and n — 1; 

2. find (r n ,r n , s n , s n ) G V£ x such that 

(r n ,s n ) G argmin -\\r<g>^y-l(r<g>s)-a(u n -i + r n <gis n ,r<Sis), 

(r,s)eV x xV t 2 

(r n ,s n )e argmin —\\r®s\\y-a(r<8>s,r n <8s n ); 

(r,s)£V x xV t 1 

(51) 

3. set u n = n n _x + r n ® s n and n — n + 1. 

The Euler equations associated to these coupled minimization problems 
read: for all (5r, 5r, 5s, 5s) G V£ x V t 2 , 

f (r n ®Sn,r n <g>8s + 5r<g>'s n )v = l{r n <g> 5s + 5r <g> s n ) - a(w n _i + r n <g> s n ,r n <g> 5s + 5r ® s n ), 
\ a(r n <g) s„,r n ® 5s + 5r ® s n ) v = a(r n ® 5s + 5r <g> s n ,r n ®s n ). 

(52) 



51 



Using Lemma [1.11 and (jlip . these definitions imply that for all n £ N*, 

n-^-n l{r® s) - a(u n -i + r n <g) s n ,r <S> s) 

\\r n ® Sn\\v = sup rp= — =n (53) 

(r,s)evix^ ||r <g» s|| v 

and _ _ 

n o n a(r <g> s, r n ® s n ) 

«Fn®Sn V= SUp . (54) 

(r,«)ev»xVt 

Indeed, for the first equality, one has to consider the symmetric coercive 
continuous bilinear form a 1 (-,-) = (-, -)y and the continuous linear form 
/ x (-) = /(•) — a(n n _! + r n (g) s„, •). For the second equality, the symmetric 
coercive continuous bilinear form to consider is a 2 (-,-) = a(-,-)y and the 
continuous linear form is l 2 (-) = a(-,r n (g) ~s n ). 
The following result holds: 

Proposition 7.1. Lei assume that the iterations of algorithm 157]) are 
well-defined for alln £ N*. Then, (u n ) n ^* converges to u in the sense of the 
infective norm: 

a(u- u n ,r®s) 
\\u-u n \\ iA = sup — >• 0. 

(r,s)ev x xv t ||r®s|| v 

Proof. Let us first remark that, using (1531 ). 

a(« — u n ,r ® s) „ 

F-\ U = SUp M-^ ~ii = F»® s n V- 

(r,3)ev*xv t F«8)s||y 

Thus, it is sufficient to prove that the sequence (||r ra (8is' Tl ||v')n€N* converges to 
as n goes to infinity. Let us first prove that this sequence is non-increasing. 
Let n > 2. From the Euler equation ( 1521) associated to the minimization 
problem defining (r n ,7 n ), it holds that 

\\r n <S> s n \\v = l(rn®s n ) - a{u n _ x + r n <g> s n ,r n ®s n ). (55) 

Besides, using (j53|) at iteration n — 1, we obtain 

n~ ~ ,, Ur <8> s) — a(w n _i,r <g> s) 
||r n _i (8) s ft _i||v = sup 

/(r n (g) g n ) - a(^ n _i, r n ® s w ) 

~ ||f n <g> sJly 
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which, using (155]) . leads to 

\\r n <8> S n \\y < ||f n _i <g> S n _i||v||f„ <8> S n ||v - «( r n ® S n ,f n S n ). 

Using the second Euler equation (1521) defining (r n , s n ), we have that a\\r n ® 
s n||y = a(r n (g) s nj r n (g> s* re ). Finally, it holds that 

||r n ® s n ||y < ||r n ® s n ||y||r n _i ® s n _i||y — oc\\r n <g> s n \\ v . (56) 

This implies that the sequence (\\r n (E) s n || v)neN* is non-increasing and thus 
converges towards a limit z > 0. Let us argue by contradiction and assume 
that 2 > 0. Dividing by \\r n (g> s re ||y equation (156|) . we obtain 

||^n <8> s n ||y ~ ~ n 

— ~ n < Fn-l S n _i v - Fn ® S n y. 
iFn ® Sn||y 

Since we have assumed that ||r n (8) S* n ||v > £ > for all uef, the series of 
general term (||r n £g> s n || y ) n6 H» converges and \\r n <E> s n ||v — ^ 0. Using (154|) . 

n— >oo 

this implies that for all (r, s) G 14 x Vt, 

a(r <g> s, f n <g> s n ) — ► 0. 

n— >oo 

Using assumption (Al) and the fact that (||r n £g> s r n ||v)neN* is bounded, we 
have 

Ww G V, a(w,r n Cg> s n ) — 0. 

n— >oo 

Since we have assumed that the operator .4. is bijective, it is surjective on V, 
and the sequence (r n (g) s n ) n , e N* weakly converges to in V. 
Using fl53D, ^ holds that for all n G N*, 

Fn®s n || y = /(f n (g)s„)-a(M n ,f n (g)s n ) = /(f n (g)? n )-a ^5Z rfe ® Sfe '^ n ® • 
Since (r n (g) s n ) ne N. weakly converges to in V, necessarily l(r n <g> s n ) — > 0. 

n— >oo 
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Besides, using (j5"4"|) . we have 



a ( <g s fc ,r n <g s n 

fc=i 



< ^ |a (r fc (g) g fc ,r ra (g s n )\ 
fc=i 

fk ® Sfe||y||r n (g) s re ||v 

fc=i 

< a I ||r fc (g) gfc||y J (n||r n <g 



2 a/2 



Since the series of general term (||r n (g s n ||y) ne ^* is convergent, the sequence 
(Efc = i Wk <g s k\\v) nm * is bounded and there exists a subsequence of (n\\r n <g 
s n \\v)n£N* which converges to 0. Thus, there exists a subsequence of (||r n (g 
S n ||v)neN* converging to and since the whole sequence converges to z, we 
have z = by uniqueness of the limit. We obtain a contradiction. □ 

Remark 7.1. Unfortunately, as announced in the beginning of this section, 
in general, the iterations of algorithm (5l\) are not well-defined in the sense 
that there may not exist a solution (r n ,r n , Snj^n) £ Kr 2 x ^t 2 °f ^ e coupled 
minimization problems. Numerically, we can observe that if we use a coupled 
fixed-point algorithm similar to the one presented for the Minimax algorithm, 
the procedure does not converge in general. Finding a suitable way to adapt 
these ideas in an implementable well-defined algorithm is work in progress. 

Remark 7.2. When the bilinear form a(-, ■) is symmetric and coercive, all 
the iterations of algorithm 157]) are well-defined. If (-,-)y is chosen to be 
equal to a(-, ■), the second equation of [5^1 implies that r n (g s n = -r n (g 
and the first equation of can be rewritten as: for all (Sr, ST) G V x x V t , 

H — (r n ®s n ,r n <g)8s+5r®s n )v = l{r n ®5s+5r® ®8s+8r®s n )v- 
a J 

This Euler equation is similar to the Euler equation of the first iteration of 
the standard greedy algorithm applied to the symmetric coercive problem: 

find u G V such that 

Vt> G V, (l + (2, u)v = - K-i, v) v . 
Thus, if we consider now the following (well-defined) algorithm 
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1. set Uq = and n = 1; 



2. find (r n , s n ) E V x x V t such that 



(r n ,s n )e argmin -\\r <g) s\\ v - l(r <g) s) - (w n -i, v) v ; 

(r,s)eV x xV t 2 



3. set u n = u n _i + r n £g> s n and n = n + 1, 

following the proof of Proposition Vfll we can prove that the sequence (ii n ) ne N* 
converges to C in the sense of the injective norm as soon as X > 1. Let us 
point out that when A = 1, i/iis algorithm is identical to the standard greedy 
algorithm applied to the symmetric coercive problem 

{find u G V such that 
W E V,{u,v) v = l(v). 



7.2 Link with a symmetric formulation 

Let us now present another approach, for which no convergence result have 
been proved so far. The idea is based on the article [7| by Cohen, Dahmen and 
Welper, where the objective was to develop stable formulations of multiscale 
convection-diffusion equations. 

The principle of the method is to reformulate the antisymmetric problem 
f lTB]) defined on the Hilbert space V, as a symmetric problem defined on the 
Hilbert space V x V. Indeed, it is proved in [7] that the unique solution of 
the problem 

find (v,v) G V x V such that 

a(w,v) = 0, VwGV, (57) 
a(v,w) — (Ryv,w) = l(w), Vu> G V, 

is (v,v) = (u, 0) where u is the unique solution of f|T6|) . 

This new problem is now symmetric. It is equivalent to the following 
problem 

find (v,v) G V x V such that 
A* \ ( v\ f 

m V x V . 



A -R v J \v J \ L 

This new formulation of the problem is symmetric, but not coercive. No 
convergence results exist for greedy methods in this framework. However, 
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the situation seems more encouraging than in the original non-symmetric 
case. It is to be noted though that the use of a simple Galerkin algorithm, 
similar to the one introduced in Section |2~3| does not work in this case either. 
More subtle algorithms need to be designed in this case as well, and this is 
currently work in progress. 
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